{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Avg_velocity_vs_depth.ipynb",
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "T_hNienaGs1z",
        "outputId": "b567ffb0-0643-4e60-dcd8-25fbe5fae0e5"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Mounted at /content/drive\n"
          ]
        }
      ],
      "source": [
        "from google.colab import drive\n",
        "drive.mount('/content/drive')"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "import pandas as pd\n",
        "from datetime import datetime\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import netCDF4 as nc\n",
        "!pip install xarray\n",
        "!pip install rioxarray\n",
        "import xarray as xr\n",
        "import rioxarray\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "fFcMzGG2W2uQ",
        "outputId": "b0b12fef-eaa6-4450-8ddb-e52c62554509"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/\n",
            "Requirement already satisfied: xarray in /usr/local/lib/python3.7/dist-packages (0.20.2)\n",
            "Requirement already satisfied: typing-extensions>=3.7 in /usr/local/lib/python3.7/dist-packages (from xarray) (4.1.1)\n",
            "Requirement already satisfied: numpy>=1.18 in /usr/local/lib/python3.7/dist-packages (from xarray) (1.21.6)\n",
            "Requirement already satisfied: importlib-metadata in /usr/local/lib/python3.7/dist-packages (from xarray) (4.12.0)\n",
            "Requirement already satisfied: pandas>=1.1 in /usr/local/lib/python3.7/dist-packages (from xarray) (1.3.5)\n",
            "Requirement already satisfied: pytz>=2017.3 in /usr/local/lib/python3.7/dist-packages (from pandas>=1.1->xarray) (2022.1)\n",
            "Requirement already satisfied: python-dateutil>=2.7.3 in /usr/local/lib/python3.7/dist-packages (from pandas>=1.1->xarray) (2.8.2)\n",
            "Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.7/dist-packages (from python-dateutil>=2.7.3->pandas>=1.1->xarray) (1.15.0)\n",
            "Requirement already satisfied: zipp>=0.5 in /usr/local/lib/python3.7/dist-packages (from importlib-metadata->xarray) (3.8.1)\n",
            "Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/\n",
            "Collecting rioxarray\n",
            "  Downloading rioxarray-0.9.1.tar.gz (47 kB)\n",
            "\u001b[K     |████████████████████████████████| 47 kB 3.3 MB/s \n",
            "\u001b[?25h  Installing build dependencies ... \u001b[?25l\u001b[?25hdone\n",
            "  Getting requirements to build wheel ... \u001b[?25l\u001b[?25hdone\n",
            "    Preparing wheel metadata ... \u001b[?25l\u001b[?25hdone\n",
            "Requirement already satisfied: xarray>=0.17 in /usr/local/lib/python3.7/dist-packages (from rioxarray) (0.20.2)\n",
            "Requirement already satisfied: packaging in /usr/local/lib/python3.7/dist-packages (from rioxarray) (21.3)\n",
            "Collecting pyproj>=2.2\n",
            "  Downloading pyproj-3.2.1-cp37-cp37m-manylinux2010_x86_64.whl (6.3 MB)\n",
            "\u001b[K     |████████████████████████████████| 6.3 MB 9.0 MB/s \n",
            "\u001b[?25hCollecting rasterio\n",
            "  Downloading rasterio-1.2.10-cp37-cp37m-manylinux1_x86_64.whl (19.3 MB)\n",
            "\u001b[K     |████████████████████████████████| 19.3 MB 1.3 MB/s \n",
            "\u001b[?25hRequirement already satisfied: certifi in /usr/local/lib/python3.7/dist-packages (from pyproj>=2.2->rioxarray) (2022.6.15)\n",
            "Requirement already satisfied: numpy>=1.18 in /usr/local/lib/python3.7/dist-packages (from xarray>=0.17->rioxarray) (1.21.6)\n",
            "Requirement already satisfied: importlib-metadata in /usr/local/lib/python3.7/dist-packages (from xarray>=0.17->rioxarray) (4.12.0)\n",
            "Requirement already satisfied: pandas>=1.1 in /usr/local/lib/python3.7/dist-packages (from xarray>=0.17->rioxarray) (1.3.5)\n",
            "Requirement already satisfied: typing-extensions>=3.7 in /usr/local/lib/python3.7/dist-packages (from xarray>=0.17->rioxarray) (4.1.1)\n",
            "Requirement already satisfied: pytz>=2017.3 in /usr/local/lib/python3.7/dist-packages (from pandas>=1.1->xarray>=0.17->rioxarray) (2022.1)\n",
            "Requirement already satisfied: python-dateutil>=2.7.3 in /usr/local/lib/python3.7/dist-packages (from pandas>=1.1->xarray>=0.17->rioxarray) (2.8.2)\n",
            "Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.7/dist-packages (from python-dateutil>=2.7.3->pandas>=1.1->xarray>=0.17->rioxarray) (1.15.0)\n",
            "Requirement already satisfied: zipp>=0.5 in /usr/local/lib/python3.7/dist-packages (from importlib-metadata->xarray>=0.17->rioxarray) (3.8.1)\n",
            "Requirement already satisfied: pyparsing!=3.0.5,>=2.0.2 in /usr/local/lib/python3.7/dist-packages (from packaging->rioxarray) (3.0.9)\n",
            "Collecting snuggs>=1.4.1\n",
            "  Downloading snuggs-1.4.7-py3-none-any.whl (5.4 kB)\n",
            "Collecting click-plugins\n",
            "  Downloading click_plugins-1.1.1-py2.py3-none-any.whl (7.5 kB)\n",
            "Collecting cligj>=0.5\n",
            "  Downloading cligj-0.7.2-py3-none-any.whl (7.1 kB)\n",
            "Requirement already satisfied: attrs in /usr/local/lib/python3.7/dist-packages (from rasterio->rioxarray) (21.4.0)\n",
            "Collecting affine\n",
            "  Downloading affine-2.3.1-py2.py3-none-any.whl (16 kB)\n",
            "Requirement already satisfied: setuptools in /usr/local/lib/python3.7/dist-packages (from rasterio->rioxarray) (57.4.0)\n",
            "Requirement already satisfied: click>=4.0 in /usr/local/lib/python3.7/dist-packages (from rasterio->rioxarray) (7.1.2)\n",
            "Building wheels for collected packages: rioxarray\n",
            "  Building wheel for rioxarray (PEP 517) ... \u001b[?25l\u001b[?25hdone\n",
            "  Created wheel for rioxarray: filename=rioxarray-0.9.1-py3-none-any.whl size=54611 sha256=7421462b46197c3e0549cea8f88256947b5adfb5e51dd8cb9cd370150d4c0f75\n",
            "  Stored in directory: /root/.cache/pip/wheels/07/da/9e/1cc57b2e7a29a206893db83e984a341e2e94378263e0798229\n",
            "Successfully built rioxarray\n",
            "Installing collected packages: snuggs, cligj, click-plugins, affine, rasterio, pyproj, rioxarray\n",
            "Successfully installed affine-2.3.1 click-plugins-1.1.1 cligj-0.7.2 pyproj-3.2.1 rasterio-1.2.10 rioxarray-0.9.1 snuggs-1.4.7\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "open=('/content/drive/MyDrive/REU_2022_copy/forcing_copy/ORA-1979-2018_all_depths_U_point_35_235.nc')\n",
        "df2=xr.open_dataset(open)\n",
        "print(df2)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "6OwPWbrhXFc9",
        "outputId": "044c42b6-3f82-4723-a3b4-2c2e5ff9ff0c"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "<xarray.Dataset>\n",
            "Dimensions:    (time: 480, LEV: 75, latitude: 1, longitude: 1)\n",
            "Coordinates:\n",
            "  * time       (time) datetime64[ns] 1979-01-15 1979-02-15 ... 2018-12-15\n",
            "  * LEV        (LEV) float64 0.5058 1.556 2.668 ... 5.698e+03 5.902e+03\n",
            "  * latitude   (latitude) float64 35.5\n",
            "  * longitude  (longitude) float64 235.0\n",
            "Data variables:\n",
            "    vozocrte   (time, LEV, latitude, longitude) float32 ...\n",
            "Attributes: (12/28)\n",
            "    cdm_data_type:              Grid\n",
            "    Conventions:                COARDS, CF-1.6, ACDD-1.3\n",
            "    creator_name:               APDRC\n",
            "    creator_type:               institution\n",
            "    creator_url:                http://apdrc.soest.hawaii.edu/thredds/dodsC/l...\n",
            "    Easternmost_Easting:        235.0\n",
            "    ...                         ...\n",
            "    standard_name_vocabulary:   CF Standard Name Table v29\n",
            "    summary:                    Reanalysis Data ORAS5 1x1 grid vozocrte opa4\n",
            "    time_coverage_end:          2018-12-15T00:00:00Z\n",
            "    time_coverage_start:        1979-01-15T00:00:00Z\n",
            "    title:                      Reanalysis Data ORAS5 1x1 grid vozocrte opa4\n",
            "    Westernmost_Easting:        235.0\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "U=True\n",
        "V=False\n",
        "salt=False\n",
        "temp=False\n",
        "ORA=True\n",
        "GODAS=False\n",
        "STD=False\n",
        "layers=True\n",
        "\n",
        "\n",
        "if U:\n",
        "  direction='u'\n",
        "  if ORA:\n",
        "    model='ORA'\n",
        "    experiment= model + '-1979-2018_all_depths_' + direction + '_235_35'\n",
        "  elif GODAS:\n",
        "    model='GODAS'\n",
        "    experiment= model + '-1980-2022_all_depths_' + direction + '_215_35'\n",
        "elif V:\n",
        "  direction='v'\n",
        "  if ORA:\n",
        "    model='ORA'\n",
        "    experiment= model + '-1979-2018_all_depths_' + direction + '_235_35'\n",
        "  elif GODAS: \n",
        "    model='GODAS'\n",
        "    experiment= model + '-1980-2022_all_depths_' + direction + '_215_35'\n",
        "elif salt:\n",
        "  direction='salinity'\n",
        "  if ORA:\n",
        "    model='ORA'\n",
        "    experiment= model + '-1979-2018_all_depths_' + direction + '_215_35'\n",
        "  elif GODAS:\n",
        "    model='GODAS'\n",
        "    experiment= model + '-1980-2022_all_depths_' + direction + '_215_35'\n",
        "elif temp:\n",
        "  direction='temp'\n",
        "  if ORA:\n",
        "    model='ORA'\n",
        "    experiment= model + '-1979-2018_all_depths_' + direction + '_215_35'\n",
        "  elif GODAS:\n",
        "    model='GODAS'\n",
        "    experiment= model + '-1980-2022_all_depths_' + direction + '_215_35'\n",
        "\n",
        "print('Running ' + experiment)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "IOkubSbgGm4L",
        "outputId": "056b77b5-025b-442a-a86f-116a0a892711"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Running ORA-1979-2018_all_depths_u_235_35\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "if U: # U \n",
        "  if ORA:\n",
        "    dir=df2.vozocrte\n",
        "  elif GODAS: # GODAS\n",
        "    dir=df2.uogrddsl \n",
        "elif V: # V\n",
        "  if ORA:\n",
        "    dir=df2.vomecrtn\n",
        "  elif GODAS: # GODAS\n",
        "    dir=df2.vogrddsl\n",
        "elif salt:\n",
        "  if ORA:\n",
        "    dir=df2.vosaline\n",
        "  elif GODAS:\n",
        "    dir=df2.a;soasdokf # don't have it yet\n",
        "elif temp:\n",
        "  if ORA:\n",
        "    dir=df2.votemper # don't have it yet\n",
        "  elif GODAS:\n",
        "    dir=df2.a;lskdfj # don't have it yet\n",
        "\n",
        "avg0=dir.mean('time')\n",
        "avg1=np.squeeze(avg0,axis=1)\n",
        "avg=np.squeeze(avg1,axis=1)\n",
        "lev=df2.LEV\n",
        "print(avg)\n",
        "print('vozocrte is', type('vozocrte'))"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "rceaXkgq1bw4",
        "outputId": "a7554dc8-3417-4601-caf8-566ab3a360be"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "<xarray.DataArray 'vozocrte' (LEV: 75)>\n",
            "array([ 0.00740551,  0.00622134,  0.005373  ,  0.0048268 ,  0.0045457 ,\n",
            "        0.00449072,  0.00463664,  0.0049676 ,  0.00549254,  0.0062449 ,\n",
            "        0.00728148,  0.00871724,  0.01072874,  0.01348455,  0.01682769,\n",
            "        0.02023705,  0.02336118,  0.02526268,  0.02604135,  0.02605538,\n",
            "        0.02520897,  0.02352811,  0.021495  ,  0.01933754,  0.01710434,\n",
            "        0.0148361 ,  0.01259707,  0.01049759,  0.00864294,  0.00707302,\n",
            "        0.0057681 ,  0.00466285,  0.00369253,  0.00278809,  0.00195748,\n",
            "        0.00117523,  0.00045154, -0.00020792, -0.00083543, -0.00140936,\n",
            "       -0.00193001, -0.0023591 , -0.0026449 , -0.00274602, -0.00265383,\n",
            "       -0.00244085, -0.00215571, -0.00187427, -0.00161672, -0.00138465,\n",
            "       -0.00122716, -0.00105003, -0.0009133 , -0.00081651, -0.00076681,\n",
            "       -0.00080029, -0.000898  , -0.00089095, -0.0008145 , -0.00068976,\n",
            "       -0.00055014, -0.00043709, -0.00035431, -0.00027486, -0.00020565,\n",
            "       -0.00018649,  0.00084785, -0.00054698,         nan,         nan,\n",
            "               nan,         nan,         nan,         nan,         nan],\n",
            "      dtype=float32)\n",
            "Coordinates:\n",
            "  * LEV        (LEV) float64 0.5058 1.556 2.668 ... 5.698e+03 5.902e+03\n",
            "    latitude   float64 35.5\n",
            "    longitude  float64 235.0\n",
            "vozocrte is <class 'str'>\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "std_positive0=(dir.std('time'))/2\n",
        "std_positive1=np.squeeze(std_positive0,axis=1)\n",
        "std_positive=avg+np.squeeze(std_positive1,axis=1)\n",
        "std_negative0=(dir.std('time'))/2\n",
        "std_negative1=np.squeeze(std_negative0,axis=1)\n",
        "std_negative=avg-np.squeeze(std_negative1,axis=1)\n",
        "# not sure if this is the correct way to show standard deviation because the orange line seems farther away from the blue line (avg) compared to the distance between the green & blue line?"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "qu1UNN61dUuv",
        "outputId": "2ffb0d87-8a07-4153-903b-fcb94d95180d"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "/usr/local/lib/python3.7/dist-packages/numpy/lib/nanfunctions.py:1671: RuntimeWarning: Degrees of freedom <= 0 for slice.\n",
            "  keepdims=keepdims)\n",
            "/usr/local/lib/python3.7/dist-packages/numpy/lib/nanfunctions.py:1671: RuntimeWarning: Degrees of freedom <= 0 for slice.\n",
            "  keepdims=keepdims)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "import matplotlib.patches as mpatches\n",
        "fig,ax=plt.subplots(figsize=(5,12))\n",
        "\n",
        "m=plt.plot(avg,lev,color=\"black\")\n",
        "\n",
        "#for legend\n",
        "#mxl=plt.axhline(y=0,color='yellow')\n",
        "\n",
        "#legend\n",
        "mxl_horiz=mpatches.Patch(color=\"#F5E663\",label='Mixed Layer')\n",
        "therm_horiz=mpatches.Patch(color='#47A8BD',label=\"Thermocline\")\n",
        "deep_horiz=mpatches.Patch(color='#FFAD69',label=\"Abyssal\")\n",
        "if STD:\n",
        "  std_col=mpatches.Patch(color='darkgrey',label=\"$\\sigma$\")\n",
        "  plt.legend(handles=[mxl_horiz, therm_horiz, deep_horiz, std_col],loc='upper left')\n",
        "else:\n",
        "  plt.legend(handles=[mxl_horiz, therm_horiz, deep_horiz],loc='upper left')\n",
        "plt.title('Average Zonal (u) Velocity Profile at (35$^{\\circ}$N, 135$^{\\circ}$W)')\n",
        "#plt.title('Average ' + direction + ' at (35N,135W) 1979-2018 ' + model) #EPGP 30-40N, 150-160W... 0-60, 120-280... 35N, 135W\n",
        "\n",
        "if salt:\n",
        "  plt.xlabel('Avg ' + direction + ' (PSU)')\n",
        "elif temp:\n",
        "  plt.xlabel('Avg ' + direction + ' (deg C)')\n",
        "elif U:\n",
        "  plt.xlabel('Avg ' + direction + ' (m/s)')\n",
        "elif V:\n",
        "  plt.xlabel('Avg ' + direction + ' (m/s)')\n",
        "\n",
        "plt.ylabel('Depth (m)')\n",
        "plt.ylim(0,1000)\n",
        "plt.gca().invert_yaxis()\n",
        "if STD:\n",
        "  s1=plt.plot(std_positive,lev, color=\"red\",linewidth=0)\n",
        "  s2=plt.plot(std_negative,lev, color=\"red\",linewidth=0)\n",
        "  plt.fill_betweenx(lev,std_positive,std_negative,color=\"darkgrey\")\n",
        "  if layers:\n",
        "    #mixed layer\n",
        "    for i in range(0,80):\n",
        "      plt.axhspan(i,i+1,facecolor='#F5E663',alpha=0.6)\n",
        "    #thermocline\n",
        "    for i in range(80,170):\n",
        "      plt.axhspan(i,i+1,facecolor='#47A8BD',alpha=0.6)\n",
        "    #abyss\n",
        "    for i in range(170,5000):\n",
        "      plt.axhspan(i,i+1, facecolor='#FFAD69',alpha=0.6)\n",
        "  plt.savefig('/content/drive/MyDrive/REU_2022_copy/plots_copy/Avg_Velo_Salt_Temp_Depth/' + experiment + '.png',format='png')\n",
        "else:\n",
        "    if layers:\n",
        "      for i in range(0,80):\n",
        "        plt.axhspan(i,i+1,facecolor='#F5E663',alpha=0.6)\n",
        "    #thermocline\n",
        "      for i in range(80,170):\n",
        "        plt.axhspan(i,i+1,facecolor='#47A8BD',alpha=0.6)\n",
        "      for i in range(170,5300):\n",
        "        plt.axhspan(i,i+1, facecolor='#FFAD69',alpha=0.6)\n",
        "plt.savefig('/content/drive/MyDrive/REU_2022_copy/plots_copy/Avg_Velo_Salt_Temp_Depth/' + experiment + '_No_STD.png',format='png')"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 732
        },
        "id": "BZdotGJ8M0Z9",
        "outputId": "55dc69b1-c523-4395-d6c9-d79808d42685"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 360x864 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWMAAALLCAYAAAAsZeP/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd3wUdf7H8dc3vUFCKAkdpIoUaVItoKAgRbAginJ2Pet59p96tju9Oz37qdjLWREFFBEQQXoH6T2hpEAghFASUr6/P3aDAZKwhN3MbvJ+Ph55ZHdmdvYzu7PvfPPdme8Yay0iIuKsIKcLEBERhbGIiF9QGIuI+AGFsYiIH1AYi4j4AYWxiIgfUBiLiPgBhbGIiB9QGFcxxpgkY8xFZcx/3hhzn4frWmiMOct71Z2w/jJrPYX1rDbGXOCFkrzKGNPKGLPcGJNtjLnn+Dq9tf3lrM3j/UBKdyqfkYAJY2PMDGNMpjEm3OlaTocx5lpjzIESfqwx5kmHa6sNXA+84+FDXgSeKWN9k40xJ8w3xgw1xqQZY0LKV+mpsdaeZa2d4X7u0wo49+MPu9+zdGPMR8aYmHKu7iHgV2ttNWvta8Xr9CUP/iCfsB8YYz4zxqQaY/YbYzYYY24uNm+GMSan2L68/rj1nWOMmWiM+WsJdewyxkQXm3azMWaGh9txlzFmsTEm1xjzUQnzS635ZHWXVLMx5lFjzE/HrWNjKdOudt8t8zNSXECEsTGmCXAuYIEhPlh/hYQCgLX2f9bamOI/wH1AOvBuRdVRij8Bk6y1hz1cfgLQxxiTWMr8j4FRxhhz3PTrgP9Za/PLV6bjBrvft05AF+Dx4xfwcJ9qDKz2cm3e8CdO3A+eB5pYa6vj+gw+Z4zpXGz+XcX26VbHre8a4Crg+OkAwcC95awzBXgO+KCU+Seruay6S6r5N6CnMSYYwBhTFwgFOh43rbl7WTj5Z+SogAhjXH+l5wMfAaOLJhpjHjbGjC2+oDHmVWPMa+7b9Ywx3xpjdhtjthpj7im2XJL78b8DB40xIcaYR4wxm93/Nq4xxgwrtnwnY8wy97xvjDFfGWOeKza/1OcqizGmI/AKcLW1NrXY9DPdf7n3uf99HVJsXpIx5gFjzO/GmCx3LRHueaVugwcGADOPq88aY5oXu/9R0XZba3OAJcDFpazve6Amrj+kRY+vAQwCPnHf9/h1O8lr0tAYM869nj3GmDeKzUsyxlxkjPkUaARMdLeEHjLGPGiM+fa453nNGPNq2S8VWGt3Aj8BbYs9z/H7VIk1G2OmA32AN9y1tDRltFhP8XUqaz8+4TUoYRUn7AfW2tXW2tyiu+6fZid7jdw+Br4EVpUw79/AA8aYOA/XVbymcdba74E9pcz3ds2LcIXv2e775wK/AuuPm7bZWpviruFkn5FjCvb7H2AT8GegM5AHJLinNwYOAdXc94OBVKA7rj80S4AngTDgDGALcLF72SRgOdAQiHRPuxKo537sCOAgUNf9+GRcf8FDgeHAEeA59+PKfK4ytisO2Aw8fNz0UPc2P+ZeX18gG2hVrPaF7lrjgbXA7WVtQ7F1JwEXlVLPbqDrcdMs0LzY/Y+Kttt9/zXgP2Vs47vAe8Xu3wYs9+R1K15rWa+J+31fAbwMRAMRQO+Stvn47Xe/vweBOPf9EGAX0LmU7Sm+roa4WrbPlrRPefA+zgBuLu29Kbp/stephBrLvQ+Uth+4p/8X1+fNAkuBmGLbsRvIAOYAF3j4uS7avnH88Vm6GZhxivnwHPBRKfNKrLm8deMK37+4b78B3Aj8/bhpHxz3mDI/I0eXO5WNduIH6I0rgGu5768r2nD3/dnA9e7b/XD9VQLoBmw7bl2PAh8W2xFuPMlzLweGAucBOwFz3PM+58lzlbJug+tfmPHF1+uedy6QBgQVm/YF8FSx2kcVm/cv4O2ytuH4D0Apy+YBrY+bdrIw/vvxO18J798+IMJ9f06xHdeT96go+Ep9TYAe7g9VSCk1FF/PCduPq3V7i/v2IGBNGduTBBxwb1Myrg97ZLF5NxZb9mTv4ww8C+NT3r/Kuw+Uth8Umxfsfk8fB0KLvY/VgHBc/7lmA808qKto+9oCWUBtvBzGpdVc3rrd+9t37tsrgBbAJcdNG33cY8r8jBT9BEI3xWhgirU2w33/c4p1Vbjvj3TfvsZ9H1yt5nrufw/3GWP24WqhJBR77PbiT2SMud64vt0uWr4tUAtXK2Ondb+yJTzWk+c63sPAWbjeuOPHMa0HbLfWFhablgzUL3Y/rdjtQ0DMSbbBE5m4ds5TUQ1XMJXIWjsbV8vjMmNMM+AcTu09KlLWa9IQSLbl74P+GBjlvj0K+PQky19mrY2z1ja21v7ZHtu3Wny/8OR99MQp7V+nuQ9AGfuBtbbA/Z42AO5wT1tgrc221uZaaz/G9Qd3oKdPZq1dBfwAPHIKNXqspJrd08tT929Ab2NMPFDbWrsRmIurLzke12v923GPKfMzUqTCvrgqD2NMJK5O9GBjTFH4hANxxpgO1toVwDfAS8aYBsAwXK0kcH0otlprW5TxFEdD0BjTGNe/1BcC86y1BcaY5bhasKlAfWOMKRacDXF1MXj6XMW36wLg/4DzrLUlvUkpQENjTFCxD3IjYMNJ1lvWNnjid6Alrr6xIoeAqGL3E4Edxe6fCXx2kvV+gqvfvxXws7U23T39VF63sl6T7UAjY0yIB4Fc0gDe3wNvGWPa4moZl9SP6qni6y/X+1gCj18nD/eBkl6D4kraD44XQun9rxbP97kif8PVjfDSKT7uVJRVM3hW9zwgFrgFV3hjrd1vjElxT0ux1m497jGefEb8vmV8GVAAtMHVQX42rg2bhevDjbV2N65/9z7EtcOudT92IZDt/kIl0hgTbIxpa4zpWspzReN6M3YDGGNuwP3FDK43oAC4y/2lzFBcLbwiHj+XcX3b+iVwn7V2WSm1LMAVgg8ZY0Ld4T3Y/biylLUNnpgEnH/ctOXANe5tuqT4fPeXhp2BqSdZ7ye4/h29BVcrtMipvEdlvSYLcf3BfMEYE22MiTDG9CqllnRcfa5HWdeXLGNxtdgXWmu3nWR7PFXe9/F4p/I6ebIPnPAaHOeY/cAYU8cYc7UxJsb93Bfj+m/0F2NMnDHmYvdrHmKMuRZXt97kU9lAa+0m4CvgmC8mjesL449Keoz7+SJwdUMEF9Vwsprd88tVt/u/oMXA/bhyqMhs97RjWsWn8Bnx7z5j9wvzUgnTr8L1b3qI+/51uHbAB49brh6uPro0XP96zafsvsO/A3tx/Vv9H1zfKN/sntcFVzAdwNUaHwc84clzHfccT7prPVDCz9vFljvL/fxZwBpgWLF5x9SOqx/rs5NtQ2nbXWxeLVyt3shi07rg+pIqG9e/71/wR1/5lcA4D9/LGe7XJby879FJXpNGuFq4e9zb/lpJ24zrO4BtuP5tfKDYMr3d78sNJ9mOsl6/kvapsmqeUdZ7c1zdHu1fHu4DJb4Gpe0HuPpyZ7qX3w+s5I8+9tq4WtDZ7vnzgX4e7hPHb29DIIdifca4wvOWUh7/FH8cJVH089TJavZC3c+7n6vTcZlkgduOW9bjz4hxP0BOkTFmAa7w/NDpWrzJGPMPYJe19hUPll0A3GRdfX4BzRjTCNeXw4nW2v1O1+O0U9kPfFhDGK4vxNpba/OcquN0nMpnRGHsIWPM+biOJ8wArgXeBs6wxY4NlsBkjAnC1YKsbq290el6pGry9z7jo4wxlxhj1htjNhljfPKt60m0wvVXeh/wV+AKBXHgM65TcffjOizybw6XI1VYQLSMjetUww24PjA7cPX1jLTWrnG0MBERLwmUlvE5wCZr7RZr7RFc30YPdbgmERGv8evjjIupz7EH0+/AdfbMUcaYW4FbAaKjwzu3alG34qqrgqy15OUVkJdfQH5+AQX5heQXFFLg/skvKKQgv+CP2+6fM1vVIzw81OnyRcpl6fKkDGttbV+sO1DC+KSstWOAMQCdOza182Y87XBFgS0r6xDrNqSwbn0KGzalkZqaSWp6Fmnp+0hL28fezIOlPjYmJoIacVHExUZTIy6auLgo4mKjiIuL5v57BlI38ZTHhBHxC+Fxo5N9te5ACeOduI5BLNLAPU28JGNPNuN/WML4H5bw+8ptpKb9cWJgaGgwdRPiSEyMo/kZCfTu2cp9P5bEhDhq16pGXJwreGOrRxIaGii7lYj/CJRPzSKghTGmKa4QvhrXOBTiBR9+MpM7//IRBQWFNDsjgYv6tKV1q3q0almXM1vXp2nj2gQHB8rXCyKBKSDC2Fqbb4y5C/gZ16mPH1hr/XFQ7oD05LNjia0eyXdf3U+3rs0wJ4wFLyK+FhBhDGCtnYTrnPlyyc8PJnVPI3LzIrxYVeXw/fiJZGUd4ogJYtnaCKKjw33SEg4PzaFuzW2EhBR4fd0igS5gwvh0pe5pRPW4+sTXqKaWXwn2Zx8mPT2L7AM55B6BajERVK8eSfVqkURGhp32+q217N2bTeoeaJhw/KBWIlJlwjg3L0JBXIbq1VzBm3sknz17ssnaf5idKZnsJJPQkGCqVYskJiac6KhwIiLDTnl8RGMM8fHV2J2h/0xESlJlwhhQEHsgPCyEenVrUK9uDfLyC8jef5j92a6fvZkHAAgKMkRFhRMd7Qrn6KhwQkODT7puvf4ipatSYSynJjQkmPj4GOLjXVeiP3Ikn4MHczl4KJeDB3PZtSuLorPpQ0KCiIwIIzKy2E9EqAJYxENVNoxzMp4De8B7KzQxRNQ64YrtxwgOrc01I6/g00/eAiA/P5/6DdtyzjmdmDj+cyZMnMzatet5+KHyXrncZcbMObz0nzeZOP5zj6Z7KiwshLCwEGrUiAZc/cCHDh/h4MFccnLyOHT4CBl7siks/GO8k4iIUHcwhxERHkpeXgGHDx/xSj+0SGVSZcPYq0Hs4fqio6NYvXothw8fJjIykqnTZlC/XuLR+UMGX8KQwZd4ty4fMsYc7aYoLjc3n/3Zh8jLK+RwjiusM91n7GVk7KNjp1tp1KAmLVok0qJZIi2aJ9LSfbthg5o6plmqpKobxg4ZMOAifpw0lSsuH8KXX37H1VcPZ9bs+QB89PEXLFmynNdf+yeXDb+O4cMGcf11I3hnzMfMmjWPzz59mylTf+Xpp/9Fbm4uZzRrwgfvvUZMTAyTf/6F+//6OFGRkfTq1e0kVRzr2ede5IcffuZwTg49unfl7bdeYsuWJEZcfROLF00HYOPGzYy85hYWL5rOkiUreODBJzhw4CA1a8Xz4fuvU7duIn0vHEqHDm2ZM2cBV189nPv/8uejz1FYaMnJzWOtzeHxh4eycXM6Gzam8tmi2WRn5xxdLjw8lObNEmjRLJEzmtahUcOaNG5Ui8YNa9GoYU2qVYv0wrsg4n8UxhVsxFXDePbvLzLo0v6sXLmaG2645mgYF/fOWy9x7vmX0rRpI15+5b/MnT2ZjIw9/OMf/2HKz2OJjo7mX/9+jZdfeZsHH7iL226/n2lTxtG8+Rlcfc3Np1TTnX++iScefwCA60f/mR9+nMLgQRcTG1ud5ctXcvbZ7fjo4y/40+iR5OXlce99j/LduE+oXbsWX339HY8/8Q/ef+81AI4cyWPhgmknPEdQkCEqMozoqHAef2TY0enWWtJ3ZbFhYxobN6excVMaGzamsnrtDib9vJwjR469vmh8jWgaNaxF40aucC663bhhLRo1qkmNuGj1U0tAUhhXsPbtzyI5aTtffDmOAQMuKnW5hIQ6PPW3R7jwomF8O/Zj4uNr8MOPU1izdgPnnncpAEfy8ujerQvr1m2kaZNGtGjhuvDttddcybvvfeJxTb/OmM2LL77BoUOH2ZuZyVlntWLwoIu58cZRfPTxF7zUrg1ffzOe+XN/Zv36TaxavZaLL7kCgIKCQhLr/nHF+KuuOrWRTY0xJCbEkZgQx3m9Wx8zr7CwkPRd+0nelsG27Rls276HpG272bZ9D+s3pjJ1+koOHTpyzGOqVYugkbsV7QrsWjSoV4O6iTWoWzeOeok1iI4+tltFxB8ojB0wePAlPPTwU0yf9j179maWutyqVWuoWTOelJQ0wNWKvOii8/n8szHHLLd8+cpy15KTk8Nddz/MwvlTadiwPk8/8y9ycnIBuHz4IJ597t/06XMunTq1P1rLWW1aM2f2TyWuLzoquty1HC8oKIi6iXHUTYyj+znNT5hvrWXP3gOusN6WQfL2PccE99x5G8naf+iEx8VWj6JuXdd66yXGUdd9KF/R7aLnDAvTx0MqjvY2B9zwp2uIi61Ou3ZtmDFzTonLLFy4lMmTf2HJoun0uXAI/ftdQPdunbn7nofZtGkLzZufwcGDB9m5M43WrVuQlLydzZu30qxZU778apzHtRQFb61a8Rw4cIBvx03k8uGDAYiIiKB/vz7cedeDvDvGdV3KVq2aszsjg3nzFtGjR1fy8vLYsGEzZ53VutTn8BVjDLVqVqNWzWp07ti0xGWysg6RkraPlNRM1zCgaftISd1Haprr9qw560lN30de3omnaNeqWY16dYta1HF/tK7r1qBeXdf9OrWr6wtH8YqqG8YmxuuHtnmqQYN63H33raXOz83N5bbb7+f9916lXr1E/v2vZ7j5lnuZNvU7Pnj/da4ddRu5ua5/z5955lFatmzG22+9xOCh1xAVGUnv3t3Jzi5526ZPn0WjJu2P3v/qi/e5+aZRtD/7PBITatOl89nHLH/NyCv4fvwk+vfrA0BYWBhff/kB9/3lMbKysskvyOeeu29zJIw9ERsbRWxsFGe2qlfqMoWFhezZe4DU1H2kpGW6fqdmkpK272iA/75yG+m7so45bA9cfeEJdWJJTIglISGOugmxJCS4hhZNrBNLYmIcCQmx1E2I0+F8UqaAuAbeqSppcPlNO86kdauSW09Supf+8yZZWft55ulHvbK+deu30rzBWq+sq6Ll5xewa/d+V1C7W9cpqZmkp2eRmp5Fevo+0tKzSgxtgOrVI4uFdKwrxBP/CO3EOq4grxkfQ1CQWtv+KDxu9BJrbRdfrLvqtozlpIZfMZotm5OYNtXzbo/KLCQk+Oip4mUpKCgkY0+2O6T3kV50hRT37/T0LJYsSyJ9VxYHDuSc8PiQkGAS6lQnMSHOHdixJNZxD+ZfFNwJ6teubPROSqnGjf3Y6RICUnBwEAl1XC3f9u0albnsgQM5pO3KIi1tH+m7skh1/05L20farix27NzLkmVb2bV7PyX9F1u7lqtfu369eOrVq0H9oi8ji27Xq0FcbJQO9wsACmMRB8XERNA8JoLmZySUuVx+fgG7M7KPhnSqu297Z6qrq2THzr0sXLyZjD3ZJzw2MjLMFdjucC5q3dev5/qpV7cGiQmxulyWw/TqiwSAkJDgo4fclSU3N4+Uoi8gUzPZmZJ59HZKaibzF24iJTXzhJNpjDEk1Kl+NKhLbGXXrUH16pFqZfuIwlikEgkPD6Vpk9o0bVL61eSLjs9OSck82rJOcYf2jpS9JCXvZu78DSVeATwmJoLGjWrRpFEtmjSuffSnaWPXfZ2uXn4KY5Eqpvjx2WX1aR8+fISUtD+CemdqJjt27CV5WwZJybuZOXvdCV9A1oyPcQf0sWHdpHFtGjesSXh4qK83L2BV2TC+5oc17MvNP/mCHooLD+HzQW1Knb9nz1769R8OQFr6LoKDg6ldqyZJydupVy+RVb+XfPKHU4oPWvT2Ox8RFRXJ9deNcLosqUCRkWE0a5pAs6Yl92cXtbCTkne7fzLY6r69YuU2Jk5adkx3iDGGenXjXOHc6MSwrl+vRpU+gabKhrE3g9iT9dWsGc/SJTMAePqZfxETE81f77+TpKRtDLns2nI/b35+PiEhvn0bb7/tTz5dvwSm4i3sLp3OOGF+YWEhKan7ioX1brYm7T7aqv7863nHHCESGhpMwwY1j7asmzauTduzGtKjWwtqxHnvNHt/VWXD2J8UFBRw621/Yd78RdSrV5fvx31CZGQkmzdv5a57HiZj9x6ioiJ55+2Xad26BTfceBcREREsX76Snj3PYe/eTCIjI1m+fCW7dmfw3ruv8umnXzF/wWLO6dqJDz94A4AvvhzHCy+8gsUycEA/Xnj+SQAm//wLjz/+dwoKCqlVK56pU449rrj4H4++Fw7lnHM6M2PGbPZlZfHumFc4t3cPCgoKePSxZ5k5cw65uUe4444bue3W0RX+Wor/CAoKokH9eBrUj6d3z1YnzM/NzWPbjj1HW9XFQ3vCD0uPOTLkrDYN6NWjJb26t6Rn95Y0alizIjelQiiM/cDGjVv436fvMOadlxkx8ia+HfcDo669ktvv+Cv/ffPftGjRjAULlnDX3Q8xbep3AOzYmcLsWZMIDg7mhhvvInPfPubM/okJEydz2bBRzJr5I2ed1Zpu3fuxfPlK6tSpzaOPPcOiBdOoUSOOSwZcyffjJ9Gr5zncdvv9zJg+gaZNG7O3jIGLiuTn5zN/3hQm/TSVZ599kSk/f8v7H/yP2NhqLJg/ldzcXM4971L697uApk0b+/rlkwAVHh7qurhAs8QS52dnH2bp8iTmzt/AnPkb+eLruYx53zW+dsMG8fR0B3PvHi1pc2b9gD9rUWHsB5o2bcTZZ7cDoHOnDiQnb+PAgQPMnbeIEVffdHS53CN/DBd5xeVDCA7+4yKggy69GGMM7dqeSUJCbdq1c/Vft2nTmqTk7SRv28H55/Widu1aAIwceTmzZs0jODiYc3v3OBqa8fFln10GMOyyS4/WmpS8DYCp035l5co1fPvtRACy9mezcdMWhbGUW7VqkZx/7pmcf+6ZgOtY61WrdzBn/gbmzt/Ab7PX8dVY11jgsdWj6NG9Bb26t6Bn95Z06dSUiIjAGgtEYewHwsP/GF83ODiYw4dzKCy0xMVVP9rPfLzo6Kjj1uHa8YKCgo5ZX1CQIT8/n9BQ732LXbT+4OBg8vNdo51Za3n1lee5uH9frz2PSHEhIcGc3aExZ3dozJ239cNaS1JyBnPmrXcH9EYmT1kBuK7X2LljE3r1aMUtN/ShSePSD/XzF4Hdrq/EqlevRtMmjflm7HjAFXYrVqwq9/rO6dqR32bNJSNjDwUFBXz11Xecd15PunfrzKzZ89i6NRnAo26KkvTv15e33/mIvLw8ADZs2MzBgycepyriLcYYmjapzaiRvXnr1RtZseB5dm5+g7Gf38tdt/fHWnj1zcl07PEYL7/+09GGg7+qsi3juPAQrx/a5m2ffvIWf77rQf7xj5fJy89jxFXD6NChbbnWVbduIv/4+xNceNGwo1/gDR0yAIC333qJK678E4WFltp1ajFl8thTXv/NN40iOXkbXbpeiMVSu1ZNxn3r+dVGRLyhVs1qDB7YicEDOwGwbfse7n3gEx554ku+/GYeb712A53O9s/RGzWEplSoQB5CUwKTtZbvJizm/oc/I31XFnfe1o+n/u9yYmIiTnldvhxCU90UIlKpGWMYPrQrKxY8z8039OH1t6bQsftj/Dh5udOlHUNhLCJVQmxsFK+/NJpfJ/8f0THhDL/6ZVau2uZ0WUcpjEWkSunZvSXvvnkzANt37nW4mj8ojEWkygkPcx3qecTLwyKcDoWxiFQ5Ye6jn3KP5DlcyR8UxiJS5RRdO/DIEf859rjKHmdsfv4XJrfky9mXhw2PwV780EmX+378JC6/YjSrV86ldesWzJg5h5f+8yYTx3/utVpOpvjwmCJVzeat6dx0+xgAEhJiHa7mD1W2ZezNID6V9X351Th69+rGl1/pissiFclay/sfz6Br7ydYszaFj8bcRr++5TuJyheqbBg74cCBA8yZs4B3x7zKV19/d3R69v5sBg0ZyZlndeeOPz9AYWEhH3z4P/5y//8dXebd9z7l/r8+zsGDBxk0ZCQdO11A+7PPPbqeRx97hrbte3F2x/N58KG/ATDxh5/p0fNiOnfpQ/+LLyc9fVfFbrCIn0hL38fwq1/hz/d+yDldmrFk7nOMvKqnX13Pr8p2Uzhh/ITJXNy/Ly1bNqNmfDxLlrgGNVm4aBmrfp9N48YNGXDpCMZ99wNXXTmU5194hX/98ylCQ0P5+OMveOu/LzL55+nUq5vIDxO+ACAraz979uzl+/GTWLNqHsYY9u3LAqB3r27MnTMZYwzvvf8p/37xDV789zOObb+IE8ZPXMyf7/uIAwdzePH5a7nztov8crhN/6uoEvvyq3GMGDEMgBFXXXa0q+Kcrh0544wmBAcHc/WIYcyZs4CYmBj6XNCbH36cwrp1G8nLy6Nduza0a9uGab/M5JFHn2HW7HnExlYnNrY6EeER3HzLvYz77geiolwXhdyxI4VLBl5Fh7PP46X/vMnqNesc23aRirZn7wFuvuNdrrrudRo1rMn8mc9w9x39/TKIQS3jCrN3bya//jqbVavWYoyhoKAAYwwDB/Y74V+lovs33TiK5//5Cq1bteBPo0cC0LJlMxYv/IVJP03jySefp2/f83ji8QeYP+9nfpn+G99+O5H//vd9pk39jnvve5T77ruDIYMvYcbMOTzzzL8qfLtFKlJBQSHTfl3FJ/+bxYQfl1JQUMhjDw7lsYeGEBrq33Hn39VVImO/ncioa6/k7bdeOjqtT98hzJ49n4WLlrF1azKNGzfk62++55abrwegW7fO7Ni+k2XLfmf50pkApKSkER8fx6hrryQuLpb3P/iMAwcOcOjQYQYO6Eevnt1o3tI1jklW1n7q16sLwCeffFnBWyxScTZsSuOT//3G/76cQ0rqPmrGx3DLDX24afQFnNWmgdPleaTKhrENj/H6oW1l+eqrcTz44N3HTBs+bBBvj/mIrl3O5u57H2Hz5q1ccH7vo1fSALjyyqEsX7GKGjXiAFi5ag0PP/w0QUGG0NBQ3nzj32RnH2DY8OvJycnFWnu0X/jJJx9ixMibqBEXS58+55KU5D/n4Yucrv37DzP2+wV88r9ZzFuwiaAgw8X92vOfF0Yx8JKzCQ/33gUVKoKG0PRzg4dew3333s6Ffc9zuhSv0BCacjoKCwv5bfZ6Pvnfb4ybsJjDh4/QqmVdrr/2XK4d0Yu6iXE+fX5fDrE/hz0AACAASURBVKFZZVvG/m7fviy69+xP+/ZnVZogFimvpOTdfPbFbD75fDbJ2zKoXj2Sa6/uxehrz6Vr5zP86hC18lIY+6m4uFjWrVngdBkijtmbeYDxE5fw5TfzmDHL9cV3n/Pb8MwTVzB0UGciIwPrgqMnU6XC2FpbKf6CBqrK2CUm3pWVdYiJk5byzXcLmTZ9Ffn5BZzRtA5/e2w4o0b2plHDmk6X6DNVJozDQ3PYuzeb+PhqCmQHWGvZuzeb8NAcp0sRP5OdfZgfJy/nm3ELmPLLSo4cyadxw1rc8+eLuXL4OXTs0KRKfGYrZRjvzonmrQ3nHDMtFEurgweJSfXumBTiuQPWsD6/FnlZ/n/ZdPGtIzmHWTN7Jsum/cSaOTPJy80ltk4CPYZfQ8d+A2jctgPGGBYACzY6XW3FqJRhXJI8DKvyA+tQF5HKJC83l7XzZrFs2k+snjWDI4cPUS2+Ft2HXEHHfgNo0r6j354dVxGqTBiLSMXLzzvC+vlzWDbtJ1b+Np3cgweJjqtBl0sG07HfAJp17EJQcLDTZfoFhbGIeFVBfj6blixgyZRJrJwxjcPZ+4msVp2z+15Mx34DadGlG8Ehip7j6RURkdNWWFhI0srlLJ3yI8t/+ZkDe/cQHh1Nu/MupGO/AbTq1pOQ0Mp1KJq3KYxFpFystezcuI6lP//IsqmTyExLJTQ8nDa9LqDTxQNp0/N8QsPDnS4zYCiMReSU7EreytIpk1g6dRK7krYQFBxC6+49ufSO+2h73oVEREc7XWJAUhiLyEllpqeybOpPLJ3yIzvWrcEYwxkdu3D+1dfToW9/YuJqOF1iwFMYi0iJDmTuZfkvk1k6ZRJbli8BoFGbdgy992E69htAXJ0EhyusXBTGInJUzoED/D5jGkun/MiGRfMoLCggsWkzBt52Dx37D6R2w8ZOl1hpKYxFqriC/DzWzZvNokkTWDVrOvlHjhBftz59Rt1I5/6XUrd5yypxOrLTFMYiVZC1lu1rV7Fo0gSWTZ3Egcy9RMfVoMfQK+l08aU0aXe2AriCKYxFqpDMtBQW/zSRRT9NYFfSFkLCwmh7bl+6DBzCmT16ExyiIQOcojAWqeRyDhxg+fSfWfzTBDYtWQjAGR060+exp+lw4SVEVavucIUCCmORSqkgP5/1C+aw6KcJrJr5C3m5udRu2JgBt91Nl0uGULN+YFyksypRGItUEtZadqxfw+JJE1g6ZRLZezOIqh5Lt8HD6TJgyNFhKcU/KYxFAty+9DQWT57I4p8mkLZlE8GhoZzV+wK6DBhCm17naUyIAKEwFglAebm5rJw5jfkTvmXjovlYa2naviNXPvw3zr7oEqJjfXuVZPE+hbFIAEndvIF548eyeNIEDu3PokbdevS/6c90GTBYJ2QEOIWxiJ/LPXSQZVN/Yt74sSSvWkFwSCjtLriQHkOvoEXXHlX66hiVicJYxA9Za0le/Tvzx49l2dRJ5B46RGLTZlx238N0GThUA/NUQgpjET9ycF8miydPZP74saRu3khYRCQd+w2g+9ArdFZcJacwFnFYYWEhmxYvYN6Esfz+61QK8vJo1KYdIx57mo4XDSQiJsbpEqUCKIxFHJK1excLf/iO+RO+Zc/O7URWq07PYSPoPuRy6rds7XR5UsEUxiIVqCA/n7Vzf2Pe+LGsmTMTW1hIi87dGHjbPbS74CLCIiKcLlEcojAWqQAHMvcy97uvmfPtF2Tt3kX1mrW48Pqb6TZ4uA5JE0BhLOJTOzesY+ZXn7L05x/IP3KEVt16ccVDT9Km1/m6XL0cQ3uDiJcVFhSw8rdf+O2rz9i8dBFhEZGcM2gY5101isQzmjtdnvgphbGIlxzan8X8Cd8y65v/kZmaQo269Rhyz4N0H3I5UdVjnS5P/JzCWOQ0pW3dzG9ffcriSRM4knOYZp26ctl9D9P23L7qihCPaU8RKYfCwkLWzZvFzC8/Zf2COYSEhdHp4kGcN2IUDVqe6XR5EoAUxiKnIOfgQRb+8B2zvv6M3duTia1dh4G330vPYVcRUyPe6fIkgCmMRTyQsWMbv339GQsmjiP34EEat23Pdbe+SIe+/TResHiFwlikDJuXLWb6Zx+wZvYMTFAwHS+6hHNHjKJJ2w5OlyaVjMJYpARbVixl8pjX2bBoPjE14ul3w230unwksbXrOF2aVFIKY5FiklYu56cxb7B+wRxi4mty2X0P03P4CMIiIp0uTSo5hbEIkLz6d34a8wbr5s0iOq4GQ+55kF6XX014ZJTTpUkVoTCWKm372tX8NOZ11syZSXRsHIPuvJ9zr7yG8Khop0uTKkZhLFXSjvVrmPzum6z6bTpR1WO59I77OPeqUUREK4TFGQpjqVJSNq5n8rtv8PuMaURWq86A2+7m/BHXawB3cZzCWKqE1M0bmPzuf1kx/WciomO4+OY7OX/k9URVq+50aSKAwlgqubQtm5j83pus+OVnwqKi6H/THVwwcrQG7hG/ozCWSik9eSs/v/smy6ZOIiwykotG38oF14wmWldVFj+lMJZKJffwISa/+yYzv/iYkNAw+l53E31G3ahL24vfUxhLpbFm7m+M/ecz7E3dSfehV3DpHfdRLb6m02WJeERhLAEve08G415+nmVTJlGnyRnc/c6nNOvYxemyRE6JwlgCVmFhIQsmjGXC6y9xJOcwl9x6FxddfwshYRpFTQKPwlgCUtrWzXz9j7+xZcUSmnXqylWPPk1C46ZOlyVSbgpjCSh5ublM+3gM0z56l/CoKEY+8XfOGTQMY4zTpYmcFoWxBIyNSxby9fN/Y/e2JDpfMpjL7ntYX9BJpaEwFr93cF8m419/kYUTx1GzfkPueP09WnXr5XRZIl6lMBa/Za1lyeSJfPfyCxzOzubC0bdw8U13aGxhqZQUxuKXMnZs45t/Ps36BXNp3LYDIx57hnrNWzpdlojPKIzF78z7/mvGvfQPgkNCueKhJ+g5bARBwcFOlyXiUwpj8RsF+Xl89/ILzP7mc1p3783Vjz9HXJ0Ep8sSqRAKY/ELB/dl8tGjf2HjkgX0HXUjg+68X61hqVIUxuK41M0beO+vd5KVsYtrn3qBrgOHOl2SSIVTGIujVv02nU+ffJDwqGjuevsTmrTt4HRJIo6olGFcxxzi3uClTpchZbDW8o8PxvP+W2PpfGZTvn/pL9SvUwDofRP/dZ8P110pw1j826HDudz4zBi+mjKfay7pyXtP3EJkhAb3kapNYSwVanvaHi77639Ytj6ZF+6+modGD9K4EiIojKUCzV2xgeEPvsKhnFwm/Od+Bp3XyemSRPxGkNMFSNXw4YSZ9Lnt71SLimD+R08riEWOo5ax+FR+fgEPvvo5r3w+mYu6teWr5+8mPjbG6bJE/I7CWHwm++BhLn/wFaYuWMW9Iy/hxfuuISREJ3KIlERhLD6Rk3uEy/76MjOXruW9J27hpssucLokEb+mMBavKygoZNQT/2X6otV8+uwdjBrY2+mSRPyevsATr7LWcuc/P+TbXxbx8l9HKYhFPKQwFq968q2xvPPtdB69YQj3XTPA6XJEAobCWLzmtS8m89z733PzZRfw9zuvcrockYCiMBav+PynOdz74qcM69OFtx69UWfViZwihbGctslzVzD6b+9wQecz+fzvd+rwNZFyUBjLaZm/ciOXP/gq7Zo3ZPx/7iciXAP+iJSHwljKbc2WHVx674vUqx3HT68/RPWYKKdLEglYCmMpl22pGfS/8wXCQkKY8uYjJNSMdbokkYBW4WFsjGlojPnVGLPGGLPaGHOve3q8MWaqMWaj+3cN93RjjHnNGLPJGPO7MUYjzDhsd+Z++t/5AgcO5fLzmw/TtH4dp0sSCXhOtIzzgb9aa9sA3YE7jTFtgEeAX6y1LYBf3PcBBgAt3D+3Am9VfMlSJPvgYQbe82+S0zL44ZUHaN+ikdMliVQKFR7G1tpUa+1S9+1sYC1QHxgKfOxe7GPgMvftocAn1mU+EGeMqVvBZQuus+tGPvYGy9Yn8fUL99C7YyunSxKpNBztMzbGNAE6AguABGttqntWGpDgvl0f2F7sYTvc045f163GmMXGmMW7M7N9VnNV9t9vpvLj7OW8fP8oBms8YhGvciyMjTExwLfAfdba/cXnWWstYE9lfdbaMdbaLtbaLrVrVPNipQKwbmsKD7zyOQN6deCuEf2dLkek0nEkjI0xobiC+H/W2nHuyelF3Q/u37vc03cCDYs9vIF7mlSQvLx8Rj3xX6Ijwnn/iVt0dp2IDzhxNIUB3gfWWmv/U2zWBGC0+/ZoYHyx6de7j6roDmQV686QCvDMu9+xZO1Wxjx+M3Vr13C6HJFKyYnxjHsB1wErjTHL3dMeA14AvjbG3AQkA0UjzUwCBgKbgEPADRVbbtU2d8UG/vHheP40+DyG9+3qdDkilVaFh7G1djZQ2v+5F5awvAXu9GlRUqLsg4e57om3aJRYi1cfuM7pckQqNV3pQ0r1l5c+Iyl1NzPHPKFTnUV8TKdDS4m+/3Ux74+fwcOjB+t4YpEKoDCWE6Rl7OOW596jU+smPHXb5U6XI1IlKIzlGNZabnrmXQ4czuHTZ+8gLFQ9WSIVQWEsx3jn21+YNGc5/7z7atqc0cDpckSqDIWxHLUhOZW/vvw5/bq11Vl2IhVMYSzAH2fZhYeF8OFTtxEUpF1DpCKpQ1AAePHTH1m0egtfv3AP9evEO12OSJWj5o+wa28Wz380gSHnd+LKft2cLkekSlIYC0+PGcehnCP8656RTpciUmUpjKu4dVtTeGfcdG4d1pdWTeo5XY5IlaUwruIeef1LoiLCeOq24U6XIlKlKYyrsFnL1jF+5hIeHj2YOvG6urOIkxTGVZS1lgde/pz6dWrwl2sGOF2OSJWnMK6ivp46n4WrN/PsHVcSFRnudDkiVZ7CuArKPZLHo298RfsWjbj+0nOdLkdE0EkfVdJ/v5nG1p27+fmNhwkO1t9jEX+gT2IVk7n/IM++9x39u7ejf4/2TpcjIm4K4yrm7+9/z77sQ/z7vmucLkVEilEYVyFbd+7i9a+m8KfB59G+RSOnyxGRYhTGVcj/vfk1wUFBPHvHFU6XIiLHURhXEYvXbOGLn+dx/6gBGpVNxA8pjKsAay0PvPI5tWtU56HrBzldjoiUQGFcBfy2dB0zl6zlyVuGUT0myulyRKQECuMq4F8f/0Cd+OrcNPQCp0sRkVIojCu5lRu3MWnOcu65+mIiI8KcLkdESqEwruRe/HQS0ZHh3HHFRU6XIiJlUBhXYtvT9vD55LncMqwP8bExTpcjImVQGFdiL3/+ExarITJFAoDCuJLK3H+QMeOmM/LiHjSqW8vpckTkJBTGldRbY6dx8HAuD+q4YpGAoDCuhHJyj/DqF5O5pGd7jUEhEiAUxpXQJz/OZtfe/TrbTiSAKIwrmYKCQl789Ee6nnUGF3Rp43Q5IuIhhXElM37mEjZuS+Oh6wdhjHG6HBHxkMK4ErHW8s+PJtKsQQLD+nR1uhwROQUK40rkt6XrWLh6Mw9cN1DXthMJMPrEViL/+vgHateozuhB5zldioicIoVxJfHHgED9NSCQSABSGFcSL346iaiIcP58ZT+nSxGRclAYVwK79mbx+eS53HzZBRoQSCRAKYwrgc8mzSG/oIDbr7jQ6VJEpJwUxgHOWssHE2bSrW0zzmxa3+lyRKScFMYBbvGaLazevIMbhpzvdCkichoUxgHuwwm/EREeytX9ezhdioicBoVxAMvJPcIXP8/l8r7nEFtNV30WCWQK4wD2/Ywl7Ms+xA1DdJKHSKBTGAewD8bPoHHdWvTR6GwiAU9hHKC2pWYwbeFqRg86l6AgvY0igU6f4gD1yY+zsNbyp8HqohCpDBTGAchay0cTZ9GnSxua1q/jdDki4gUK4wA0a9k6Nu9I1xd3IpWIwjgAfTB+JtWiI7i87zlOlyIiXqIwDjDZBw/zzbSFjOjXnajIcKfLEREvURgHmG+mLeBQTq5OfxapZBTGAebDCb/RqnFderRv4XQpIuJFCuMAsiE5ldnL13PDkPN15WeRSkZhHEA+mvgbQUGG6y7t7XQpIuJlCuMAUVBQyCc/zuKSHh2oV7uG0+WIiJcpjAPELwtXsXNXpo4tFqmkFMYBYuwvC4mJimDQuR2dLkVEfEBhHAAKCgr5fsYSLu19NhHhYU6XIyI+oDAOAHNWrGd35n6G9+nqdCki4iMK4wAwbvpiwsNCGdCrg9OliIiPKIz9nLWWcdMX0b97O6pFRzpdjoj4iMLYzy1Zu5Xt6XsY3qeL06WIiA8pjP3cuOmLCA4OYvB5nZwuRUR8SGHsx6y1fDt9ERd0PpOacdWcLkdEfEhh7MfWbt3JhuRUHUUhUgUojP3Yd78uBuAy9ReLVHoKYz82bvoierRvobEoRKoAhbGfSkrZzdJ1SeqiEKkiFMZ+qqiLYpi6KESqBIWxnxo3fREdWjaiWcMEp0sRkQqgMPZDaRn7mLNig7ooRKoQhbEfGj9zCdZahvdVGItUFQpjPzRu+iJaNErkrGYNnC5FRCqIwtjPZO4/yPRFaxjep6suOipShSiM/cwPs5aSX1CgLgqRKkZh7GfGTV9Mg4R4urRp6nQpIlKBFMZ+5ODhHCbPW8GwC7oQFKS3RqQq0Sfej0xftIac3Dwuu0AneohUNQpjPzJtwSoiw8PodXZLp0sRkQqmMPYj0xau4rxOrQkPC3W6FBGpYApjP5GyO5M1W3Zy0TltnS5FRBygMPYT0xetBuCibmc5XImIOEFh7CemLVhFrbhqtG/RyOlSRMQBCmM/YK1l2sJV9O3aRoe0iVRR+uT7gfVJqezclan+YpEqTGHsB6YtXAXARd0UxiJVlcLYD0xbsIoz6tehaf06TpciIg5RGDssP7+AXxevUatYpIpTGDts8dot7D94WP3FIlWcwthh0xaswhhDny5tnC5FRBykMHbYtIWr6diqMbVqVHO6FBFxkMLYQQcP5zB3xQb1F4uIwthJs5atJy+/QP3FIqIwdtIvC1cTHhZK77NbOV2KiDhMYeygaQtX0bN9CyIjwpwuRUQcpjB2yO7M/Sxfn6z+YhEBFMaOOTpkpvqLRQSFsWOmLVhFbEwUnc/UVaBFRGHsCGstUxe4hswMDtZbICIKY0ds2bGL5NQMdVGIyFEKYwdoyEwROZ7C2AEzl6ylXu0atGiU6HQpIuInFMYOmL9qEz3bt8AY43QpIuInFMYVbNfeLLbu3E23ts2dLkVE/IjCuIItWLUZgG7tmjlciYj4E4VxBVuwahPBwUF0bq3ji0XkDwrjCjZ/5SbaN29EVGS406WIiB9RGFeggoJCFq7eTPd26i8WkWMpjCvQuqQUsg/m0K2t+otF5FgK4wq0YNUmALWMReQECuMKNH/lJuKqRelkDxE5gcK4Ai1YtZlubZsTFKSXXUSOpVSoIAcO5bBq83b1F4tIiRTGFWTxmi0UFlqdeSciJVIYV5CjZ96pZSwiJVAYV5D5KzfRvGECNeOqOV2KiPghhXEFsNYyf+UmHdImIqVyLIyNMcHGmGXGmB/c95saYxYYYzYZY74yxoS5p4e7729yz2/iVM3ltT1tD2l79qm/WERK5WTL+F5gbbH7/wRettY2BzKBm9zTbwIy3dNfdi8XUIr6i9UyFpHSOBLGxpgGwKXAe+77BugLjHUv8jFwmfv2UPd93PMvNAE2Kvv8lRsJDwulfYtGTpciIn7KqZbxK8BDQKH7fk1gn7U2331/B1Dffbs+sB3APT/LvfwxjDG3GmMWG2MW787M9mXtp2zBqs10PrMJYaEhTpciIn6qwsPYGDMI2GWtXeLN9Vprx1hru1hru9Su4T9HLOTl5bNk3Vb1F4tImZxoqvUChhhjBgIRQHXgVSDOGBPibv02AHa6l98JNAR2GGNCgFhgT8WXXT6/b9pOTm6eji8WkTJVeMvYWvuotbaBtbYJcDUw3Vp7LfArcIV7sdHAePftCe77uOdPt9baCiz5tMxfuRHQl3ciUjZ/Os74YeB+Y8wmXH3C77unvw/UdE+/H3jEofrKZcGqzSTUjKVRYi2nSxERP+boN0rW2hnADPftLcA5JSyTA1xZoYV50fyVm+jetjkBdgCIiFQwf2oZVzp79mWzcVua+otF5KQUxj60cLVO9hARzyiMfWjJ2q0AdD6zqcOViIi/Uxj70IoN2zijfh2qx0Q5XYqI+DmFsQ+t2LiNs1s1droMEQkACmMfOXAoh03b0+mg8ShExAMKYx9ZuWk71lo6tFQYi8jJKYx9ZMWGZAB1U4iIRxTGPrJiwzbiqkXpzDsR8YjC2EeWb0imfYtGOvNORDyiMPaBwsJCVm7aztkt1UUhIp5RGPvA5h27OHg4V1/eiYjHFMY+sHy968u7DmoZi4iHFMY+sGJjMsHBQZx1Rv2TLywigsLYJ1Zs2EbrJvWICA9zuhQRCRAKYx9YsWGbzrwTkVOiMPayvVkH2J6+Ryd7iMgpURh7WdGZd2oZi8ipUBh72YqN2wAdSSEip0Zh7GXL1yeTWDOOhJqxTpciIgFEYexlKzZu08keInLKFMZedCQvnzVbdiqMReSUKYy9aF1SCkfy8jUmhYicMoWxFx09kkItYxE5RQpjL1qxYRvhYaG0bFTX6VJEJMAojL1oxcZttGvekJCQYKdLEZEAozD2Emsty9cn62QPESkXhbGXpGbsI2NftvqLRaRcFMZeoguQisjpUBh7SdGA8u2bq2UsIqdOYewlKzZuo0m92sRWi3K6FBEJQApjL1mxYZtO9hCRclMYe8GRvHw2bk+jbbMGTpciIgFKYewFW3bsoqCgkFZNdLKHiJSPwtgL1ienANCqscJYRMpHYewF65NTAXQatIiUm8dhbIyJNsboPN8SrE9OJaFmrI6kEJFyKzWMjTFBxphrjDE/GmN2AeuAVGPMGmPMv40xzSuuTP+2PilVXRQiclrKahn/CjQDHgUSrbUNrbV1gN7AfOCfxphRFVCj31ufrDAWkdMTUsa8i6y1ecdPtNbuBb4FvjXGhPqssgCxN+sAGfuyFcYiclpKDePiQWyMqQE0LL68tXZpSWFd1WzY5vryTmEsIqejrJYxAMaYZ4E/AZsB655sgb6+KytwrE9yh3GTeg5XIiKB7KRhDFwFNLPWHvF1MYFofXIqIcHBNK1X2+lSRCSAeXJo2yogzteFBKr1yak0a1BHV/cQkdPiScv4eWCZMWYVkFs00Vo7xGdVBZD1yak6DVpETpsnYfwx8E9gJVDo23ICS0FBIZu2pzOwVwenSxGRAOdJGB+y1r7m80oCUHJqBrlH8mjVWF/eicjp8SSMZxljngcmcGw3xVKfVRUgNECQiHiLJ2Hc0f27e7FpOrSNPwYIUp+xiJyuk4axtbZPRRQSiDYkp1GjejS14qo5XYqIBLiyBgoaZYwpa34zY0xv35QVGIrGpDDGOF2KiAS4slrGNXEd0rYEWALsBiKA5sD5QAbwiM8r9GPrk1O4sGtbp8sQkUqgrLEpXjXGvIGrb7gX0B44DKwFrrPWbquYEv3TgUM57NyVqf5iEfGKMvuMrbUFwFT3jxSzIVkDBImI9+iyS+W0XmEsIl6kMC6n9cmpGGNo3jDB6VJEpBJQGJfT+uRUmtSrRUR4mNOliEgl4Ml4xuHA5UATjh1c/hnfleX/dN07EfEmT1rG44GhQD5wsNhPlWWtZcO2VI1JISJe48np0A2stZf4vJIAkrI7k4OHc9UyFhGv8aRlPNcY087nlQSQoksttWyc6HAlIlJZlNoyNsasxDUgUAhwgzFmC65R2wxgrbXtK6ZE/6PD2kTE28rqphhUYVUEmPXJqURHhlO/TrzTpYhIJVHW6dDJAMaYT6211xWfZ4z5FLiuxAdWAeuTU2jZSAMEiYj3eNJnfFbxO8aYYKCzb8oJDLrunYh4W1lDZD5qjMkG2htj9htjst33d+E63K1KysvLJzk1gxYN9eWdiHhPqWFsrX3eWlsN+Le1trq1tpr7p6a19tEKrNGvbE/fS2GhpWn92k6XIiKViCfHGT9mjBkO9MZ1dMUsa+33vi3LfyWl7AagSd1aDlciIpWJJ33GbwK3AyuBVcDtxpg3fVqVH0tOywCgST21jEXEezxpGfcFzrTWWgBjzMfAap9W5ceSUnYTFGRooMPaRMSLPGkZbwIaFbvf0D2tSkpKyaB+7XhCQz35OyYi4hlPEqUasNYYsxBXn/E5wGJjzAQAa+0QH9bnd5JSd9OknvqLRcS7PAnjJ31eRQBJTs3gvE6tnS5DRCqZk4axtXamMaYx0MJaO80YEwmEWGuzfV+ef8nPL2DHrr00TlTLWES866R9xsaYW4CxwDvuSQ2AKnlo245deykoKNSRFCLidZ58gXcn0AvYD2Ct3QjU8WVR/uroMcbqMxYRL/MkjHOttUeK7hhjQnB9kVfl/HHCh1rGIuJdnoTxTGPMY0CkMaYf8A0w0bdl+afktAyMMTRMrOl0KSJSyXgSxo8Au3GdgXcbMAl43JdF+auklAzq1Y4jTMcYi4iXeXI0RaEx5nvge2vt7gqoyW8lpexWF4WI+ERZQ2gaY8xTxpgMYD2w3hiz2xhTZY87TkrN0JEUIuITZXVT/AXXURRdrbXx1tp4oBvQyxjzlwqpzo/k5xewI30vjTVam4j4QFlhfB0w0lq7tWiCtXYLMAq43teF+ZuU3ZnkFxRo6EwR8YmywjjUWptx/ER3v3Go70ryT0mpRccYq5tCRLyvrDA+Us55lVJSisYxFhHfKetoig7GmP0lTDdAhI/q8VtFJ3w0TNA4xiLifaWGsbU2uCIL8XfJaRnUrRVHRHiY06WISCXkyUkfgqubQl0UIuIrCmMPuU740JEUIuIbCmMPFBQUsi1tj44xFhGfURh7IDXDfYyxuilExEcUxh44elibxqUQER9RGHvgRp4+XAAAFj5JREFUjxM+1E0hIr6hMPZA0THGjXTtOxHxEYWxB5JTM0ioGUtkhI4xFhHfUBh7ICk1Q4e1iYhPKYw9kJSyW0dSiIhPKYxPorCwkOTUDB1JISI+pTA+idSMfeTlF+iEDxHxKYXxSSSnFg2dqTAWEd9RGJ9E0WFt6jMWEV9SGJ9Ekrtl3FjHGIuIDymMT2J72h5qxsYQFRnudCkiUokpjE8ibU8WibXinC5DRCo5hfFJpO/NIrFmrNNliEglpzA+ibQ9+0isqZaxiPiWwrgM1lrS9+wnIb6606WISCWnMC7DgUM5HMrJVZ+xiPicwrgMaXuyAEiIV5+xiPiWwrgM6e4w1hd4IuJrCuMypO3ZB6BuChHxOYVxGdL3qptCRCqGwrgMaRlZBAUZasVVc7oUEankFMZlSNuzjzo1YgkO1sskIr6llClD+t79JNTUMcYi4nsK4zLo7DsRqSgK4zKkZWTpyzsRqRAK41JYa12DBNVSGIuI7ymMS7Ev+xBH8vLVTSEiFUJhXIr0o6dC6ws8EfE9hXEpjp59p5axiFQAhXEpigYJUp+xiFQER8LYGBNnjBlrjFlnjFlrjOlhjIk3xkw1xmx0/67hXtYYY14zxmwyxvxujOlUETWma8Q2EalATrWMXwUmW2tbAx2AtcAjwC/W2hbAL+77AAOAFu6fW4G3KqLAtD37CA0Jpkb16Ip4OhGp4io8jI0xscB5wPsA1toj1tp9wFDgY/diHwOXuW8PBT6xLvOBOGNMXV/Xmb53P3XiqxMUpJ4cEfE9J5KmKbAb+NAYs8wY854xJhpIsNamupdJAxLct+sD24s9fod72jGMMbcaYxYbYxbvzsw+7SJ19p2IVCQnwjgE6AS8Za3tCBzkjy4JAKy1FrCnslJr7RhrbRdrbZfaNU5/lLW0DF0VWkQqjhNhvAPYYa1d4L4/Flc4pxd1P7h/73LP3wk0LPb4Bu5pPpW+N4sEhbGIVJAKD2NrbRqw3RjTyj3pQmANMAEY7Z42Ghjvvj0BuN59VEV3IKtYd4ZPFBYWuk6FVjeFiFSQEIee9274//buPsayur7j+Ps7986duwvyXI2AEYw0xmpMZYMkYmoFdbUPmGhVSisxpMS0TRONbWmtwVjT1LSNaWNiQ6OtWis+1Up8pvjQSqIIFgQawVVQQRFlYYfZnX2Y2W//OL/ZOb3MLrMzc8/vsvN+JZO5c+65Z345ZD9853vO73f4UEQMgO8Dr6f5H8NHI+Jy4AfAq8u+nwVeDuwA9pR9x+rBXXMsLh509p2kzlQJ48y8Bdi2wlsXrrBvAn8w9kG1LD+I1MpYUje8b2sFzr6T1DXDeAXOvpPUNcN4BS4SJKlrhvEK7n9wF8OZaU44fkvtoUjaJAzjFfx0Z/O4pYioPRRJm4RhvIL7f/6ws+8kdcowXsFPd846+05SpwzjFbhIkKSuGcYjFhYW+dlDj9imkNQpw3jEzx9+hMz0HmNJnTKMRxy6x/g02xSSumMYjzg0Fdo2haQOGcYjHnx4DoDTTlr/AvWStFqG8Yhdc3sAOPH4rZVHImkzMYxHzO6eB+CE45wKLak7hvGI2d3z9Hs9hjPTtYciaRMxjEfM7p7nxOO3uC6FpE4ZxiN2ze1xtTZJnTOMR8zunrdfLKlzhvGI2TnDWFL3DOMRu+bmva1NUucM4xG2KSTVYBiPMIwl1WAYj5jdPe/dFJI6Zxi37Nt/gH37D9gzltQ5w7jFqdCSajGMW2bnDGNJdRjGLUuV8Yn2jCV1zDBuWVo+08pYUtcM45ZDPWMv4EnqmGHc4gU8SbUYxi3LT/kwjCV1yzBu8W4KSbUYxi2zu+eZ7veYGfiUD0ndMoxbltal8CkfkrpmGLe4fKakWgzjFhcJklSLYdzi8pmSajGMW3zkkqRaDOOWXXN77BlLqsIwbrFNIakWw7jFMJZUi2Fc7N23n/0HFpwKLakKw7hYXrHNMJbUPcO4cMU2STUZxoWLBEmqyTAulpfP9NY2Sd0zjAvbFJJqMowLL+BJqskwLqyMJdVkGBe7ygU8e8aSajCMi9m5eQbTfZ/yIakKw7hwKrSkmgzjwhXbJNVkGBdWxpJqMowLH7kkqSbDuLAyllSTYVzM7dnL8Vtmag9D0iZlGBcHFhYZTPdrD0PSJmUYFwcWFun3erWHIWmTMoyLAwuLTPcNY0l1GMbFwuIifcNYUiWGcWFlLKkmw7hYWFyk3/N0SKrD9Cmayti7KSTVYRgDmcni4kHbFJKqMYyBhYVFANsUkqoxfWhaFIBtCknVGMY0F+/AylhSPaYP7crYnrGkOgxjYGHxIGCbQlI9hjFwYGEBsE0hqR7TB9sUkuozjFluU7g2haRaDGOW2xRWxpJqMYzxAp6k+gxjlnvGXsCTVIvpg20KSfUZxsDCQrmA52OXJFViGGNlLKk+w5j2BTzDWFIdhjHtC3iGsaQ6DGNsU0iqzzDGGXiS6jOMcW0KSfUZxiwvLm8YS6rFMMYLeJLqM4xZfiCplbGkWgxjrIwl1WcY4wU8SfUZxngBT1J9hjG2KSTVZxhjZSypPsOY5cq45+LykioxfWjCuN/rERG1hyJpkzKMadoUtigk1WQYUyrjvqdCUj0mEM2qbT4ZWlJNhjHNesY+GVpSTSYQTZvCnrGkmgxjmqdDO+FDUk2GMU2bwspYUk2GMUsX8AxjSfUYxixP+pCkWgxjnPQhqT7DGCtjSfUZxnhrm6T6DGNJmgCGMdDvTbF48GDtYUjaxAxjmid8LCwaxpLqMYyBfn/q0NM+JKkGw5hSGS8YxpLqMYxpesa2KSTVZBgDvSnbFJLqMoyBfr/n3RSSqjKMWeoZG8aS6jGMWeoZ26aQVI9hjPcZS6rPMMbKWFJ9hjHQM4wlVWYY07QpFhez9jAkbWJVwjgi3hgRd0TE7RHx4YgYRsTZEfGNiNgRER+JiEHZd6b8vKO8f9ZGj8c2haTaOg/jiDgD+CNgW2Y+C+gBrwXeCbwrM58OPARcXj5yOfBQ2f6ust+Gai7gGcaS6qnVpugDWyKiD2wFfgK8CPh4ef/9wCvK64vLz5T3L4yI2NDB9Kc4eDA56MQPSZV0HsaZeR/wt8APaUJ4F3Az8HBmLpTd7gXOKK/PAH5UPrtQ9j919LgRcUVE3BQRN/3soUeOakxLj1xa9PY2SZXUaFOcTFPtng2cDhwHbF/vcTPz6szclpnbfuHkJxzVZ3tTzWnwXmNJtdRoU1wE3J2ZP8vMA8C/A88HTiptC4AzgfvK6/uApwCU908EHtzIAR2qjG1TSKqkRhj/EDg/IraW3u+FwP8CXwZeVfa5DPhUeX1t+Zny/pcyc0PvQ+v3lypjL+JJqqNGz/gbNBfivgXcVsZwNfCnwJsiYgdNT/i95SPvBU4t298EXLnRY1qqjF0sSFIt/cfeZeNl5lXAVSObvw+ct8K+e4HfGud4+j0rY0l1OQOP9gU8w1hSHYYxzeLy4N0UkuoxjFluU3ifsaRaDGNaF/BsU0iqxDCmHcZWxpLqMIzxbgpJ9RnGNIvLg5WxpHoMY1woSFJ9hjG2KSTVZxjj3RSS6jOMaS0U5NoUkioxjLEyllSfYYyLy0uqzzBmeW0KF5eXVIthTOtuigXbFJLqMIxxOrSk+gxjvM9YUn2GMd5NIak+wxjXppBUn2GMi8tLqs8wxjaFpPoMY7ybQlJ9hjGttSmsjCVVYhjTmg7tQkGSKjGMsWcsqT7DmNbdFK5NIakSw5j2fcZWxpLqMIyBqakppqbCnrGkagzjot/rWRlLqsYwLnpTU95nLKkaw7jo96esjCVVYxgX0/0+B1xcXlIlhnExHEyzb/+B2sOQtEkZxsVwMM1ew1hSJYZxMZyZZn7f/trDkLRJGcbFcDDN3n1WxpLqMIyL4YxtCkn1GMbFcDCwMpZUjWFcbLEyllSRYVzYppBUk2Fc2KaQVJNhXDT3GXtrm6Q6DONiOOOtbZLqMYyL4WCaecNYUiWGcdFcwNtPZtYeiqRNyDAuhoNpDh5MFly5TVIFhnExHAwAvL1NUhWGcbFlOA0YxpLqMIyLQ5WxF/EkVWAYF8OBlbGkegzjYjhTwtg1jSVVYBgXVsaSajKMi6XK2Kd9SKrBMC4OVcZewJNUgWFcDGe8z1hSPYZxYWUsqSbDuPACnqSaDONiy9A2haR6DONiuU3h3RSSumcYF7YpJNVkGBczXsCTVJFhXPR6U0z3e1bGkqowjFuGM9POwJNUhWHcMhwMbFNIqsIwbhkOpm1TSKrCMG5pHkpqGEvqnmHcsmXGNoWkOgzjlqZN4QU8Sd0zjFuGM9NWxpKqMIxbvIAnqRbDuMULeJJqMYxbhoNp5vfaM5bUPcO4ZTgYWBlLqsIwbrFNIakWw7hlOPBuCkl1GMYtVsaSajGMW7bMDNi3/wCZWXsokjYZw7hl6Wkf+6yOJXXMMG4ZzvjoJUl1GMYtQx+9JKkSw7hlKYx92oekrhnGLcOZAWCbQlL3DOMW2xSSajGMW7yAJ6kWw7hluTK2ZyypW4ZxyxZ7xpIqMYxbDrUp7BlL6phh3HKoTWFlLKljhnGLF/Ak1WIYtwwHTc/Yp31I6pph3OIMPEm1GMYtW4elMjaMJXXMMG6Znu7T7/XYY5tCUscM4xFbhwPDWFLnDOMRW4cz7Nm7r/YwJG0yhvEIK2NJNRjGI5owtjKW1C3DeETTprAyltQtw3iElbGkGgzjEVbGkmowjEdYGUuqwTAeYWUsqQbDeIS3tkmqwTAeYZtCUg2G8YilNkVm1h6KpE3EMB6xdTggM9nnAvOSOmQYj9g6nAGwbyypU4bxiKU1je0bS+qSYTzCylhSDYbxCCtjSTUYxiOsjCXVYBiP2LplqTI2jCV1xzAesVwZ26aQ1B3DeMRyz9jKWFJ3DOMRXsCTVINhPMILeJJqMIxHWBlLqsEwHrFlxp6xpO4ZxiOmpqYYzkxbGUvqlGG8Ap/2IalrhvEKfNqHpK4ZxitoKmPbFJK6M7Ywjoj3RcQDEXF7a9spEXFdRHy3fD+5bI+I+IeI2BER346I57Y+c1nZ/7sRcdm4xttmZSypa+OsjP8F2D6y7Urg+sw8B7i+/AzwMuCc8nUF8B5owhu4CngecB5w1VKAj5PPwZPUtbGFcWb+F7BzZPPFwPvL6/cDr2ht/0A2vg6cFBFPBl4KXJeZOzPzIeA6Hh3wG84LeJK61u/49z0pM39SXt8PPKm8PgP4UWu/e8u2w21/lIi4gqaqBpiLcy+9c72DjXMvXe8h2k4Dfr6RBxwzxztejne8xjXep47hmED3YXxIZmZEbNgjmDPzauDqjTreRouImzJzW+1xrJbjHS/HO16Pt/FC93dT/LS0HyjfHyjb7wOe0trvzLLtcNsl6ZjSdRhfCyzdEXEZ8KnW9teVuyrOB3aVdsYXgJdExMnlwt1LyjZJOqaMrU0RER8GXgicFhH30twV8dfARyPicuAHwKvL7p8FXg7sAPYArwfIzJ0R8ZfAN8t+b8/M0YuCjxcT20I5DMc7Xo53vB5v4yUyN6xtK0laI2fgSdIEMIwlaQIYxmsUEdsj4s4yhfvKFd6fiYiPlPe/ERFntd77s7L9zoh46WqPOYHjvScibouIWyLipkkYb0ScGhFfjoi5iHj3yGfOLePdUabfx4SP9yvlmLeUrydOwHhfHBE3l/N4c0S8qPWZSTy/Rxrv2M7vmmSmX0f5BfSA7wFPAwbArcAzR/b5feAfy+vXAh8pr59Z9p8Bzi7H6a3mmJM03vLePcBpE3Z+jwMuAN4AvHvkMzcC5wMBfA542YSP9yvAtgk7v78MnF5ePwu4b8LP75HGO5bzu9YvK+O1OQ/YkZnfz8z9wDU0U7rb2lO/Pw5cWCqFi4FrMnNfZt5NcwfJeas85iSNd5zWPN7M3J2ZXwP2tncu97WfkJlfz+Zf4gdYno4/ceMds/WM938y88dl+x3AllKVTur5XXG8GzSuDWUYr81qpmkf2iczF4BdwKlH+Oyqp35PyHgBEvhi+fPvCjbOesZ7pGPe+xjHXKtxjHfJP5c/od+6gX/2b9R4Xwl8KzP38fg4v+3xLhnH+V2TatOhdUy4IDPvK7226yLiO9ksEKWNcWk5v08APgH8Lk3FWV1E/BLwTpqJWBPvMOOdqPNrZbw2q5mmfWifiOgDJwIPHuGz45z6PY7xkplL3x8APsnGtS/WM94jHfPMxzjmWo1jvO3z+wjwb0zI+Y2IM2n+e78uM7/X2n8iz+9hxjvO87smhvHafBM4JyLOjogBzQWDa0f2aU/9fhXwpdJLuxZ4bemznU2zhvONqzzmxIw3Io4rFQURcRxNxXE7G2M9411RNtPrZyPi/PLn6OtYno4/ceONiH5EnFZeTwO/zgSc34g4CfgMcGVm3rC086Se38ONd8znd21qX0F8vH7RTN++i+Yq71vKtrcDv1leD4GP0VzwuhF4Wuuzbymfu5PWFeeVjjmp46W5sn1r+bpjwsZ7D81a2nM0/cVnlu3baP7BfQ94N2UG6iSOl+Yui5uBb5fz+/eUu1hqjhf4C2A3cEvr64mTen4PN95xn9+1fDkdWpImgG0KSZoAhrEkTQDDWJImgGEsSRPAMJakCWAYayJFxCsiIiPiGR3/3oiIL0XECUf5uemI+NYR3r8mIs5Z/wh1rDKMNakuAb5Wvnfp5cCtmTl7lJ+7ALjhCO+/B/iTNY9KxzzDWBMnIo6nCbfLaWZbLa1n+7HWPi+MiE+X15dHxF0RcWNE/FOMrAtc9nlbRLy59fPt0VqzueVSysyxiDgrIr4TEf9Sjv+hiLgoIm6IiO9GRHv67Hbgc2Vm4mci4tbyO15T3v9v4KIyVVd6FMNYk+hi4POZeRfwYEScC/wn8Lwy9RrgNcA1EXE68FaadXSfD6y3rfF8mplZS54O/F057jOA36b5H8WbgT9v7ferNOvjbgd+nJnPycxnAZ8HyMyDNLPDnrPO8ekYZRhrEl1Cs2Yt5fsl2SyL+HngN0p1+Ws0Fex5wFczc2dmHqCZErsep2SzcMySuzPzthKmdwDXZzNt9TbgLICIOAPYmZl7yvYXR8Q7I+IFmbmrdawHgNPXOT4do/yTSRMlIk4BXgQ8OyKS5ikPGRF/TBPMf0izjsNNmfnIUSxBu8D/Lz6Gh9svIqZK+AK017492Pr5IMv/frYDXwDIzLsi4rk0ved3RMT1mfn21u+cX+2AtblYGWvSvAr4YGY+NTPPysynAHcDLwC+CjwX+D2WK+dvAr8SESeXivmVhznuPeWzlLA8+zD73UmzCNLR2E7zmCFK22RPZv4r8DdLv7P4RWqvDKaJZRhr0lxCs/Zs2ydoWhWLwKeBl5XvZLMm7V/RrNR1A03o7uLRPgGcEhF30FTXdx3m938GeOFqBxsRPeDpmfmdsunZNEuM3gJcBbyj7PckYD4z71/tsbW5uGqbHvci4vjMnCuV8SeB92XmaKCv9lhPBj6QmS9e5f4XAL+TmW94jP3eCMxm5nvXMi4d++wZ61jwtoi4iKYn+0XgP9Z6oMz8Sbk97oTV3GuczcNEv7aKQz8MfHCt49Kxz8pYkiaAPWNJmgCGsSRNAMNYkiaAYSxJE8AwlqQJ8H9bbQThxYobNgAAAABJRU5ErkJggg==\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Calculation of Terminal Velocity\n",
        "\n",
        "See ipad for calculation... Assume unit cube with density of 1.40g/cm<sup>3</sup>"
      ],
      "metadata": {
        "id": "NRoVhKcuIgnv"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Calculation of Amount of Time to Reach Terminal Velocity\n",
        "\n",
        "I just followwed [this youtube video](https://www.youtube.com/watch?v=E1ZEZ_UVCXc) for the most part\n",
        "\n",
        "## This ended up being for turbulent flow... see part_sink.ipynb for updated calculations"
      ],
      "metadata": {
        "id": "Hy_yj0KylYzm"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Assume that the density of water is constant at 1.0 g/cm^3= 10^6 g/m^3.\n",
        "\n",
        "We want to know when velocity reaches terminal velocity (maximum speed). We have calculated that, for a PET plastic unit cube (0.01m x 0.01m x 0.01m), the terminal velocity is 0.209272 m/s. PET has a density of 1.40*10^6 g/m^3, and since ρ =mass/volume, 1.40*10<sup>6</sup> g/m<sup>3</sup> = mass/10<sup>-6</sup>m<sup>3</sup>. So mass=1.40g. \n",
        "\n",
        "Let C<sub>d</sub> = drag coefficient for a square = 1.28\n",
        "<br> \n",
        "Let A=area of cube (0.01m*0.01m)= 0.0001m<sup>2</sup>\n",
        "\n",
        "\n",
        "We account for 3 different forces acting on the plastic:\n",
        "\n",
        "\n",
        "1.   Gravity (downward force): F<sub>G</sub>= mass*gravity\n",
        "2.   Drag (upwawrd force): F<sub>D</sub> = 1/2 * ρ<sub>cube</sub> v<sup>2</sup> * C<sub>d</sub> * A\n",
        "3.   Buoyant Force (upward force): F<sub>B</sub> = ρ<sub>water</sub> * (volume of displaced water) * gravity\n",
        "\n",
        "Then, we create a while loop to determine velocity of particle over time. We can solve for Net Force so that we can find acceleration since F<sub>net</sub>=mass*acceleration. Then solve for velocity and then for position. When velocity=terminal velocity, the while loop should stop and output the position (height) of the particle at that time. This gives how much far the particle must fall in order to reach terminal velocity.\n",
        "\n",
        "**I don't know how to check the accuracy of this, as it seems as though the code is only running through one iteration of the loop. I also don't know any physics, so there may be some problems with how I accounted for the forces acting on the plastic. It should take longer than 0.01 seconds to reach terminal velocity??**\n",
        "\n",
        "\n"
      ],
      "metadata": {
        "id": "Cfh_eeQJ-ykj"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# We don't need to use this\n",
        "v_terminal=-0.209272\n",
        "#v_terminal=-0.209272\n",
        "m=1.40 #grams\n",
        "g=-9.81 # or 9.81? m/s^2\n",
        "l=0.01 #m\n",
        "rho_water=10.26*(10**5) # g/m^3\n",
        "#rho_cube=1.04*rho_water\n",
        "C_d=1.28 # based on seawater & cube\n",
        "\n",
        "y=0.0\n",
        "v=0.0\n",
        "t=0.0\n",
        "dt=0.1\n",
        "\n",
        "#F=(0.5*rho_cube*(v**2)*1.28*0.0001)-(m*g)+(rho_water*0.000001*g)\n",
        "#print(F)\n",
        "while v>v_terminal:  \n",
        "  F=((l**3)*rho_water*g)+(0.5*C_d*(v**2)*rho_water*(l**2))-(m*g)\n",
        "  #F=(0.5*rho_water*(v**2)*1.28*0.0001)-(m*g)+(rho_water*0.000001*g)\n",
        "  a=F/m\n",
        "  v=v+a*dt\n",
        "  print(F)\n",
        "  y=y+v*dt\n",
        "  t=t+dt\n",
        "  print(\"depth v_terminal is reached =\", -y, \"m\")\n",
        "  print(\"time to reach terminal =\", t, \"seconds\")\n",
        "  # does this even make sense? it feels like it is only completing one iteration of the loop"
      ],
      "metadata": {
        "id": "wKYZsFF-legi"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "For a 0.01m x 0.01m cube"
      ],
      "metadata": {
        "id": "6KQ_-nmrZRbJ"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "g=9.81 #m/s^2\n",
        "print('g is of type', type(g))\n",
        "l=0.01 #m\n",
        "print('l is of type', type(l))\n",
        "rho_water=1025000.0 #g/m^3\n",
        "print('rho_water is of type', type(rho_water))\n",
        "C_d=1.28\n",
        "print('C_d is of type', type(C_d))\n",
        "m=1.40 #g\n",
        "print('m is of type', type(m))\n",
        "\n",
        "v_terminal=np.sqrt(2.0*[(1.4*9.8)-((0.01**3.0)*1025000*9.8)]/(1.28*1025000*(0.01**2.0)))\n",
        "#print(v_terminal)\n",
        "# according to wolfram alpha, v_terminal=0.209272"
      ],
      "metadata": {
        "id": "hGaDShcYoOlW"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        ""
      ],
      "metadata": {
        "id": "6NgMeSUTY_Br"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}